% RSTS_MHRUN.M
%	samples from the posterior density using Metropolis-Hastings
%  Proposal Density : Mixture of normal and t-distribution  
%
%
% Sungbae An
% created: 06/29/2004
% updated: 08/17/2005
% modified: 08/23/2006 by Taeyoung Doh  
% updated:  12/13/2007 by Taeyoung Doh 
%--------------------------------------------------------------------
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;


% MHRUN
%	odd integer  (e.g., 1) for linear posterior draws
%	even integer (e.g., 2) for quadratic posterior draws
mhrun = 1;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------


append  = 1;		% append mode
diagcov = 0;		% use diagonal covariance using previous run stat


% scaling factors, number of blocks
     
    cc     = 0.25;   % m212s5t1p1    * (0.5 and cc0: 2, without trend) 
	cc0    = 1;		% m212s5t1p1 in searching for starting point
	nblock = 250;	% number of simulation blocks
	nsimul = 10000;	% number of simulation draws in each block


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;
global psel 


 lfile = [path_.post 'pm.mat']; 
%lfile = [path_.para 'stat_mhrun1.mat'];
load(lfile); 
%pmod = drawmean'; 
path_.mhrun = [path_.para 'mhrunpmchain2' num2str(mhrun)  '/']; 
% lfile2 = [path_.post 'paranew.mat']; 
% load(lfile2); 
% pmod = paranew; 
chkdir(path_.mhrun);  
 
%pmod = meand;
% global variables for parameter specification

% parameter names
pnames = feval(fnames_.pnames,msel);


%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);


%--------------------------------------------------------------------
% data and parameters
%--------------------------------------------------------------------

% load data
loaddata;

% starting (or true) paremeters
para = pmod.*prior_.maskinv + prior_.fix.*prior_.mask;

% posterior mode and covariance matrix
%	pmod    : posterior mode
%	pmax    : posterior evaluated at the mode
%	plik    : likelihood evaluated at the mode
%   ppri    : prior evaluated at the mode
%   pmod_std: std of posterior mode estimates
%	sigmult : square root of inverse Hessian

%% use diagonal covariance matrix based on previous run
%sigmult = -sigmult; 
%% use diagonal covariance matrix based on previous run
sigmult = -sigmult; 
covd    = sigmult*sigmult'; 
%covd = drawsig; 
corrd   = zeros(size(covd,1),size(covd,2)); 
  for i=1:size(covd,1);
      for j=1:size(covd,2);
          corrd(i,j) = covd(i,j)/(sqrt(covd(i,i))*sqrt(covd(j,j))); 
      end 
  end 

para = para-sqrt(diag(covd));  
sbar        =0.9; 
df          =5; 


nobs  = size(data,1);	% number of observations
npara = length(para);	% number of parameters



%--------------------------------------------------------------------
% starting point
%--------------------------------------------------------------------
if append		% append mode

	% load para_app.mat
	lfile = [path_.mhrun 'para_app.mat'];

	if exist(lfile,'file')
		load(lfile);

		paranew = paraapp';
		postnew = postapp;
		likenew = likeapp;

		if draws==nsimul
			iblock0 = blocks+1;
			iblockn = nblock;
			j0      = 1;
		elseif draws<nsimul
			iblock0 = blocks;
			iblockn = nblock;
			j0      = draws+1;

			lfile = [ path_.mhrun 'mhdraw' num2str(iblock0,'%04.0f') '.mat'];
			load(lfile);
		end

	
		% print cotinuing point
		fmt='%15.8f ';
		disp(' ');
		disp(' ');
		disp(sprintf('Metropolis-Hastings continues at ... '));

	else
		append  = 0;
	end
end

if ~append		% new MH chain mode

	iblock0 = 1;
	iblockn = nblock;
	j0      = 1;

	%--------------------------------------------------------------------
	% find starting seed
	%--------------------------------------------------------------------
	start=0;

	% set random state for filtering or load previous one
	
	

	while ~start
		% set time
		t1 = clock;
        
			% evaluate the log-likelihood by Kalman filter
		   
		      llik = tloglik4(para,data);
       
	    % print results on screen
	    fmt='%17.6f ';
	    disp(' ');
	    disp('------------------------------------------------------------');
		disp(sprintf('Model Specification:  %s',mpspec));
		disp(sprintf(['MHRUN              :  %d'],mhrun));
	    disp('------------------------------------------------------------');
		disp(' ');
	    disp('likelihood comparison');
	    disp('------------------------------------------------------------');
	    disp(sprintf(['log-likelihood by Kalman filter   :' fmt],llik));
		disp(sprintf(['Time elapsed                      :' fmt 'seconds'],etime(clock,t1)));
		disp(' ');
	    disp(' ');

		start=1; 
	end

	% posterior value with starting seed
	ppri0 = logprior(para);
	plik0 = llik;

	%--------------------------------------------------------------------
	% look for a starting point
	%--------------------------------------------------------------------
	search=1;
	if search==1

		% print head materials
		fmt='%15.8f ';
		disp(' ');
		disp(' ');
		disp(sprintf('----------------------------------------------------------------------'));
		disp('Parameter  Fixed          Values  ');
		for i=1:length(para)
			disp(sprintf(['%s      %d '  repmat(fmt,1,1) ],pnames(i,:),prior_.mask(i),para(i)));
		end
		disp(' ');
		disp(sprintf(['Log Prior      : ' fmt],ppri0));
		disp(sprintf(['Log Likelihood : ' fmt],plik0));
		disp(sprintf(['Log Posterior  : ' fmt],ppri0+plik0));
		disp(sprintf('----------------------------------------------------------------------'));
		disp(' ');
		disp(' ');
		disp(sprintf('looking for a starting point ...'));
		disp(sprintf('----------------------------------------------------------------------'));


		% loop until finding a appropriate starting point
		% log posterior eavaluated at a starting point is greater than POSTCUT
		postcut = -1000;

		% initialize loop
		i = 1; postmax=-1e6; postnew=postmax;
		while (postnew < postcut)
		 	paranew = para + cc0*mvnrnd(zeros(npara,1),covd)';
            paranew = paranew.*prior_.maskinv + prior_.fix.*prior_.mask;

			
			     [postnew,likenew]= rsfnmh(paranew,data);
            

			if postnew > postmax
				postmax = postnew;
				disp(sprintf(['iter : %7d      log posterior: ' fmt],i,postnew));
			end
			i = i+1;
		end

	else
		paranew = para;

		
	        [postnew,likenew]= rsfnmh(paranew,data);   
        
	end

	% save starting point
	sfile = [path_.mhrun 'paranew.mat'];
	save(sfile,'paranew','likenew','postnew');

	% print starting point
	fmt='%15.8f ';
	disp(' ');
	disp(' ');
	disp(sprintf('Metropolis-Hastings starts at ... '));

end  % append

% print starting point
disp(sprintf('----------------------------------------------------------------------'));
disp('Parameter  Fixed          Values  ');
for i=1:length(para)
	disp(sprintf(['%s      %d '  repmat(fmt,1,1) ],pnames(i,:),prior_.mask(i),paranew(i)));
end
disp(' ');
disp(sprintf(['Log Prior      : ' fmt],postnew-likenew));
disp(sprintf(['Log Likelihood : ' fmt],likenew));
disp(sprintf(['Log Posterior  : ' fmt],postnew));
disp(' ');
disp(sprintf(['Scaling Factor : ' fmt],cc));
disp(sprintf('----------------------------------------------------------------------'));
disp(' ');


%--------------------------------------------------------------------
% Metropolis-Hastings simulation at a starting point, PARANEW
%--------------------------------------------------------------------

% set clock
t0 = now;

% update
likeold = likenew;
postold = postnew;
paraold = paranew;

% block loop
for iblock=iblock0:iblockn;

	% set time
	t1 = now;

	% initialize block output
	if j0==1
		likesim = zeros(nsimul,1);
		postsim = zeros(nsimul,1);
		parasim = zeros(nsimul,npara);
		rej     = zeros(nsimul,1);
	end

	% within-block loop
	for j=j0:nsimul

		% propose new parameters
		% paranew ~ N( paraold, cc^2*inv(-pmhess) )
        s = rand; 
        if s<=sbar 
        paranew = paraold + cc*mvnrnd(zeros(npara,1),covd)';
        else 
        paranew = paraold + cc*diag(diag(covd))*mvtrnd(corrd,df)'; 
        end  
		paranew = paranew.*prior_.maskinv + prior_.fix.*prior_.mask;

		    if psel == 1
			modsize = 1000;		% for reporting purpose
            else 
            modsize = 500;
            end 
            
	        [postnew,likenew]=rsfnmh(paranew,data);
       

		%----------------------------------------------------------------------
		% accept/reject decition
		% case 1: (postnew > postold)
		%	take the candidate as new draw
		% case 2: (postnew < postold)
		%   take the candidate as new draw with probability (posterior odd)
		%----------------------------------------------------------------------
		r = min([1 ; exp(postnew-postold)]);

		if (rand<r)	% accept

			% save draws in output matrices
			likesim(j)   = likenew;
			postsim(j)   = postnew;
			parasim(j,:) = paranew';

			% update
			likeold     = likenew;
			postold     = postnew;
			paraold     = paranew;

		else			% reject, discard new parameters and keep old parameters

			% save draws in output matrices
			likesim(j)   = likeold;
			postsim(j)   = postold;
			parasim(j,:) = paraold';
			rej(j)       = 1;

		end


		% print log
		if mod(j,modsize)==0
			disp(' ');
			disp(' ');
			disp(' ');
			diary on;
			disp(' ');
			disp(sprintf('block: %2.0f / %2.0f        simulation steps: %5.0f / %5.0f',iblock,iblockn,j,nsimul));
			disp(sprintf('-----------------------------------------------------'));
			disp(sprintf(['r                      : ' fmt],r));
			disp(sprintf(['Acceptance Ratio       : ' fmt],1-sum(rej)/j))
			disp(' ');
			disp(         'Parameters             : ');
			for i=1:length(para)
				disp(sprintf(['     %s           '  fmt ],pnames(i,:),parasim(j,i)));
			end
			disp(' ');
			disp(sprintf(['Likelihood             : ' fmt],sum(likesim)/j));
			disp(sprintf(['Posterior              : ' fmt],sum(postsim)/j));
			disp(' ');

			disp(        ['Elapsed Time (Block)   :  ' datestr(now-t1,13)]);
			disp(		 ['Elapsed Time (Total)   :  ' datestr(now-t0,'dd') ' days,  ' datestr(now-t0,13) ]);
			disp(sprintf('-----------------------------------------------------'));
			diary off;

			%save results in the meanwhile
			sname = [ path_.mhrun 'mhdraw' num2str(iblock,'%04.0f') '.mat'];
			save(sname, 'likesim', 'postsim', 'rej', 'parasim');

			% save for append mode
			blocks  = iblock;
			draws   = j;
			paraapp = parasim(j,:);
			likeapp = likesim(j);
			postapp = postsim(j);

			sfile = [path_.mhrun 'para_app.mat'];
			save(sfile,'paraapp','likeapp','postapp','blocks','draws');

		end

	end	% within-block loop

	%save results for each block
	sname = [ path_.mhrun 'mhdraw'  num2str(iblock,'%04.0f') '.mat'];
	save(sname, 'likesim', 'postsim', 'rej', 'parasim');

	% save for append mode
	blocks  = iblock;
	draws   = nsimul;
	paraapp = parasim(end,:);
	likeapp = likesim(end);
	postapp = postsim(end);

	sfile = [path_.mhrun 'para_app.mat'];
	save(sfile,'paraapp','likeapp','postapp','blocks','draws');

	j0 = 1;
        
end		% block loop



%--------------------------------------------------------------------
% save and print
%--------------------------------------------------------------------

% calculate elapsed time and print into log
%diary on;
disp(' ');
disp('-------------------------------------------------------');
disp(['Start       :  ' datestr(t0)]);
disp(['End         :  ' datestr(now)]);
disp(['Elapsed Time:  ' datestr(now-t0,'dd') ' days,  ' datestr(now-t0,13) ]);
disp('-------------------------------------------------------');
%diary off;

